Phase transition and selection in a four-species cyclic Lotka-Volterra model 
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We study a four species ecological system with cyclic dominance whose individuals are distributed on a 
square lattice. Randomly chosen individuals migrate to one of the neighboring sites if it is empty or invade 
this site if occupied by their prey. The cyclic dominance maintains the coexistence of all the four species if 
the concentration of vacant sites is lower than a threshold value. Above the treshold, a symmetry breaking 
ordering occurs via growing domains containing only two neutral species inside. These two neutral species can 
protect each other from the external invaders (predators) and extend their common territory. According to our 
Monte Carlo simulations the observed phase transition is equivalent to those found in spreading models with 
two equivalent absorbing states although the present model has continuous sets of absorbing states with different 
portions of the two neutral species. The selection mechanism yielding symmetric phases is related to the domain 
growth process whith wide boundaries where the four species coexist. 

PACS numbers: 05.50.+q, 87.23.Cc 



Multispecies ecological models with spatial extension ex- 
hibit a large variety of possible stationary states as well as 
phase transitions when tuning the model parameters. In the 
original Lotka-Volterra models |2[ as well as in the gener- 
alized versions the spatial distribution of species is neglected 
(see 0101 for reviews). Now we report a phenomenon under- 
lying the role of spatial effects in the biological evolution. 

In the simplest spatial Lotka-Volterra models the individ- 
uals of competitive species are residing on the sites of a lat- 
tice and the system evolution is governed by invasions along 
the nearest neighbor links. In many cases the species form 
domains whith growing sizes and sooner or later only one 
species will survive. Significantly different behavior is found 
if the species dominate cyclically each other, i.e., the corre- 
sponding food web is characterized by a directed ring graph 
Frachebourg and Rrapivsky |6] have shown that fix- 
ation occurs if the number of species N s exceeds a threshold 
value Nf(d) depending on the spatial dimension d. In this 
case the species form a frozen domain structure |6]. Con- 
versely [N s < Nf(d)], the moving invasion fronts maintain 
a self-organizing polydomain structure. These patterns are 
widely studied for N„ = 3 |5|,|8|,|9] because it can provide 
a stability against some external invaders for the spatial mod- 
els (13, [TTl E3l - Sato et al. have shown that, if only one 
of the invasion rates differs from unity for even N s , then only 
the species with odd (even) labels survive . Very recently, the 
species biodiversity were studied by similar models in bacte- 
rial (3> phytoplankton 1 15] systems. 

In the above lattice models each site is occupied by an in- 
dividual of the competitive species. Now we will consider 
a diluted version of these models for N s — 4. Namely, the 
sites may be empty and the individuals are allowed to jump 
to these empty sites. These elementary events can result in 
the formation of "defensive alliances" consisting of two neu- 
tral species. These two-species mixed states can preserve their 
territory from the external invaders belonging to the remaining 



two species. Thus, beside the above mentioned four-species 
state, this model has two sets of "defensive alliances" whose 
confrontations will determine the final stationary state. When 
increasing the concentration of vacant sites this system un- 
dergoes a phase transition from the symmetric four-species 
state to one of the symmetric defensive alliances. This tran- 
sition will be interpreted by considering the average displace- 
ment (and velocity) of boundary separating two competitive 
domains. 

In the present model the site i of a square lattice can be 
empty (si — 0) or single occupied by one the four species (i.e., 
Si = 1, 2, 3, and 4) dominating cyclically each other (1 beats 
2 beats 3 beats 4 beats 1). The time evolution is controlled 
by subsequent jumps or invasions at randomly chosen nearest 
neighbor sites i and j. The individual will jump to the empty 
site, i.e., the value of s, and Sj are exchanged (s* <-> Sj) if 
S{ = and Sj > or Sj > and Sj = 0. Invasion occurs if 
predator and prey meet. For example, both the (1,2) and (2,1) 
pairs transform into (1,1) [the further elementary invasions are 
given by cyclic permutation of the species labels]. Nothing 
happens if — Sj as well as for neutral pairs, i.e., the pairs 
(1,3), (3,1), (2,4), and (4,2) remain unchanged. The system is 
started from a random initial state. After some transient time 
the system reaches a stationary state we study. 

Notice that the above elemantary rules leave the number of 
vacant sites unchanged and their distribution becomes uncor- 
rected after a suitable relaxation time. The states containing 
only one species are considered as absorbing states because 
the above rule does not create new species. Besides this, the 
mixed states containing only two neutral species (1+3 or 2+4) 
are also absorbing states and will be denoted as and £>24- 
In these stationary states the ratio of the two species remains 
constant. In the presence of vacant sites the migration elimi- 
nates the spatial correlations. For small sizes this system can 
easily reach one of these absorbing states and afterwards it 
stays there forever. 
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FIG. 1: Spatial distribution of the four species on the square lattice 
if pa — 0. The grayscale of the four species are indicated above the 
snaphot. 

Our Monte Carlo (MC) simulations are performed on a 
square box with periodic boundary conditions. In order to 
avoid the above mentioned small size effect the linear size is 
varied from L = 400 to 2000. The systematic simulations are 
started from a random initial state for different concentration 
of vacant sites (po). Within a time unit (MCS) each pair has 
a chance once on the average to modify the state at one of the 
corresponding sites. During the simulations we have recorded 
the concentration of species and the pair configuration prob- 
abilities on the nearest neighbor sites. Averaging over a suit- 
able sampling time interval we have determined the average 
species concentrations (p a , a = 1, 2, 3, and 4). Furthermore, 
we have deduced two quantities P pp and P n describing the 
probability of finding predator-prey and neutral pairs on two 
nearest neighbor sites. Evidently, P pp measures the invasion 
activity that vanishes in the absorbing states. 

The visualization of species distribution shows a self- 
organizing polydomain structure in the absence of vacant sites 
(for a typical snapshot see Fig. 0. The species occur cycli- 
cally at each sites, however, the short range interaction can 
not synchronize these accidental events. In the pattern evo- 
lution one can easily recognize the traveling invasion fronts 
that play crucial role in the maintenance of this polydomain 
structure |5, 9]. Similar spatiotemporal patterns can be ob- 
served for low concentration of vacant sites. Henceforth this 
spatiotemporal pattern is called C state. 

In the stationary C state p\ = p2 = pj, = Pa = (1 — Po)/4 
due to the cyclic symmetry. Strikingly different behavior oc- 
curs if po > p cr = 0.0623(1). When using lighter (darker) 
grayscales for the species 1 and 3 (2 and 4) two types of 
growing domains (namely D13 and D24) can be distinguished 
as shown in Fig. |2] These growing domains are separated 
by wide regions of C states. The growth process is similar 




FIG. 2: Typical domain structure at time t = 3000 MCS (Monte 
Carlo steps per sites) if initially (t — 0) the spatial distribution was 
random for po = 0.1. The white boxes refer to empty sites while the 
grayscale of species as in Fig.Q 



to those observed in systems with two equivalent absorbing 
states Finally the present system develops 

into one of the symmetric two-species absorbing states -D13 
or D24 where px = p 3 = (1 - p )/2 and p 2 = Pi — 0, or 
P2 = Pi = (1 — Po)/2 and p\ = p% = 0. The time of transi- 
tion toward one of these states depends on p and L. Both the 
species 1 and 3 benefit their spatially mixed coexistence be- 
cause they protect each other from the external invasions. For 
example, the species 2 can invade the sites occupied originally 
by species 3, however, the neighboring species 1 strikes back 
and eliminates the invaders 2. At the same time, species 3 
protects 1 against 4. This is the reason why this association is 
called defensive alliance. Due to the cyclic symmetry species 
2 and 4 can form a similar defensive alliance. 

The formation of defensive alliances was already observed 
in some other multispecies Lotka-Volterra model where the 
cyclic invasion itself has provided the protection mechanism 
I1 1 L 1 1 2fl . In the present model, however, the protection is due 
to the mixing of neutral species via the jumps to empty sites. 

Figure [5] demonstrates that the probability of neutral pairs 
(P n ) increases with the concentration of vacant sites in the 
stationary states. Above the mentioned threshold value (po > 
p cr ) this quantity tends to the uncorrected value, P n = 
(1 — po) 2 /2, characteristic to the symmetric defensive al- 
liance state. Simultaneously, the invasion activity (or P pp ) 
decreases and drops suddenly to zero at po = p cr . The dis- 
continuous transition is accompanied with enhanced fluctu- 
ations (in all the quantities we studied) and a critical slow- 
ing down. To avoid the undesired effects of fluctuation en- 
hancements, the above MC data have been obtained on large 
system (L — 2000) with long relaxation and sampling times 
(t r > 10 4 MCS and t s > 10 5 MCS) in the close vicinity of 
the transition point. 
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FIG. 3: Monte Carlo results for the probability of finding predator- 
prey (closed diamonds) and neutral pairs (open squares) on two near- 
est neighbor sites in the stationary states. 
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FIG. 4: Average displacement of boundaries between the states C 
and D13 as a function of time for po = 0.056, 0.060, 0.064, and 
0.070 (from top to bottom). 

0.03 r 



The universal behavior of the nonequilibrium transitions 
into absorbing states have extensively been studied for sev- 
eral decades (for a review see the Refs. 12(1 |2JJ] ). The dy- 
namical systems with two equivalent absorbing states repre- 
sent a curious universality class (named after the voter model) 
II (a ll 7l ll 8l l22l 12311 whose general features are consistent to 
those described above. 

To have a deeper insight into the dynamics of the present 
model we now will study the displacement of interfaces sepa- 
rating the C state and one of the stationary defensive alliances 
(Z?i3 or Z?24)- For the preparation of such an artificial do- 
main structure the whole area (torus) is divided into parallel 
strips with width of 500 lattice units. The MC simulation is 
started from a random initial state (as above) for L — 4000 
and after a suitable relaxation time t r uncorrected D13 state 
is created in every second strip. More precisely, species 1 and 
3 are substituted randomly for the occupied sites located in- 
side the corresponding strips. First we consider the results 
obtained for symmetric distribution, i.e., when inside the de- 
fensive alliances (Z) 13 ) px = p 3 = (1 — p )/2- The expansion 
(or shrinking) of the C domains can be monitored by eval- 
uating the quantity = pi(t) - p 2 (t) + p 3 (t) - pi{t). 
Notice that $ vanishes (($) = 0) for the state C whereas 
= ±(1 — po) in the absorbing states. The average dis- 
placement (measured in lattice unit) of the parallel interfaces 
are derived straightforwardly from the variation of 

Figure|4]displays the typical time dependences of the aver- 
age displacement d(t) of the boundaries separating the C and 
D13 states. The increase of d(t) corresponds to the expan- 
sion of C domains. The MC data are obtained by averaging 
over 20 runs performed for L = 4000 and t r = 3000 MCS 
if po < p cr . Above the critical point (po > p cr ) we have 
to use significantly shorter relaxation times (t r = 200 MCS) 
to avoid the difficulties caused by the appearance of the D13 
and D24 nucleons. It is remarkable that at the beginning d(t) 
decreases suddenly. After a suitable transient time, however, 
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FIG. 5: Average velocity of invasion front between the states C and 
D13 as a function of the concentration of vacant cites. The arrow 
indicates the critical point derived from the investigation of the sta- 
tionary states. The statistical error is comparable to the symbol size. 



the variation of d(t) becomes linear and the fitted slope can be 
interpreted as the average velocity v of the invasion front. 

Figure [5] clarifies that the C state invades the territories 
of defensive alliances for low concentration of vacant sites. 
The average invasion velocity decreases monotonously with 
po and becomes zero at po — p C r- In agreement with the ex- 
pectation, the area of the C domains shrinks for po > p C r- 

The above simulations were repeated by choosing asym- 
metric compositions (e.g., p\ > p 3 ) within the D13 state. It is 
found, that the asymmetry has influenced only the short time 
behavior. For example, if p\ ^S> P3 then the C state can invade 
fast (v ~ 1) those neighboring patches occupied by only the 
species 1 in the D13 domains. Consequently, in this case one 
can observe a sudden increase (instead of decrease as plotted 
in Fig. [4} in d(t). In the subsequent linear region, however the 
average velocity v becomes independent of (pi — p 3 ) within 
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the statistical error. The visualization of the species distribu- 
tion has indicated that the boundary between the C and Z) 13 
domains fluctuates very intensively. In fact, it is not a well de- 
fined boundary because the sites occupied by species 1 and 3 
can belong to both phases within a boundary layer. Within this 
boundary layer the cyclic invasions sustain the equivalence 
between pi and p% on average for long times. This boundary 
layer can be considered as a symmetric species reservoir that 
drives an equalization between the different species concen- 
tration in the asymmetric Z? 13 phase via diffusion. 

This scenario is checked by considering the evolution from 
such an initial state where the parallel strips (created as above) 
are filled alternately by the symmetric D24 and asymmetric 
-D13 states. The simulations (for po > p cr ) have confirmed 
that the variation of d(t) is similar to a random walk mean- 
time the difference p% — p% tends to zero for long times. This 
is the reason why we have always found symmetric D13 or 
D24 states after the domain coarsening process for sufficiently 
large system sizes. At the same time this phenomenon can be 
interpreted as a selection mechanism favorizing the symmetric 
defensive alliances. 

Here it is worth mentioning that the traditional mean-field 
and pair approximations (for details see 1 20, 24]) are capable 
to reproduce the existence of the above mentioned absorbing 
states. However, these techniques are not capable to describe 
the observed phase transition. We think that this failure is due 
to the very complex mechanisms consisting of many elemen- 
tary steps within a local cycle. 

In summary, our work shows that a slight migration in 
the lattice Lotka-Volterra models may significantly affect the 
species biodiversity. Dilution and migration are attributes usu- 
ally found in ecosystems whose description should include 
these features. The present model exemplifies that the mi- 
gration of species supports the formation defensive alliances 



in the multispecies ecological systems. Furthermore, the con- 
frontation between the different associations of species plays 
crucial role in the selection of the survival population struc- 
ture. In this case a phase transition occurs when the rate of 
migration is increased by allowing more and more vacant sites 
(and jumps) on the lattice. 

The above described features are observed for many other 
systems. Preliminary results indicate clearly that similar be- 
havior occurs if the mixing is provided by the site exchange 
for neutral pairs without introducing vacant sites. Further- 
more, quantitatively similar behavior is found for the contin- 
uous version of the present model, i.e., when the individuals 
move freely on a planar surface and they create an offspring 
if they eat a prey caught within a short distance. In fact, this 
former finding inspired us to introduce a simpler model for the 
more rigorous analysis. 

We think that the mixing of neutral species can results in 
other defensive alliances in multispecies systems. For ex- 
ample, two equivalent alliances are expected to emerge in 
the A s -species model with a circular food web for even N s . 
Evidently, such alliances can occur for more complicated 
food webs when the species have several preys and preda- 
tors 00. In these situations the competition between the 
possible (defensive) alliances will affect the evolution of the 
ecological system including the food web itself 1251 12611 . 
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